Fitting models in Pumas

2021-04-20
using Random
using Pumas
using PlottingUtilities
using PumasPlots

Fitting a PK model

In this tutorial we will go through the steps required to fit a model in Pumas.jl.

The dosage regimen

We start by simulating a population from a two compartment model with oral absorption, and then we show how to fit and do some model validation using the fit output.

The dosage regimen is an dose of 100 into Depot at time 0, followed by two additional (addl=2) doses every fourth hour

repeated_dose_regimen = DosageRegimen(100, time = 0, ii = 4, addl = 2)

1 rows × 10 columns

timecmtamtevidiiaddlratedurationssroute
Float64Int64Float64Int8Float64Int64Float64Float64Int8Route…
10.01100.014.020.00.00NullRoute

As ususal, let's define a function to choose body weight randomly per subject

choose_covariates() = (Wt = rand(55:80),)
choose_covariates (generic function with 1 method)

and generate a population of subjects with a random weight generated from the covariate function above

pop =  map(i -> Subject(id = i,
                events = repeated_dose_regimen,
                observations = (dv = nothing,),
                covariates = choose_covariates()),
                1:24)
Population
  Subjects: 24
  Covariates: Wt
  Observables: dv

We now have 24 subjects equipped with a weight and a dosage regimen.

The PK model of drug concentration and elimination

To simulate a data set and attempt to estimate the data generating parameters, we have to set up the actual pharmacokinetics (PK) model and simulate the data. We use the closed form model called Depots1Central1Periph1 which is a two compartment model with first order absorption. This requires CL, Vc, Ka, Vp, and Q to be defined in the @pre-block, since they define the rates of transfer between (and out of) the compartments

mymodel = @model begin
  @param   begin
    cl  RealDomain(lower = 0.0, init = 1.0)
    tv  RealDomain(lower = 0.0, init = 10.0)
    ka  RealDomain(lower = 0.0, init = 1.0)
    q   RealDomain(lower = 0.0, init = 0.5)
    Ω   PDiagDomain(init = [0.9,0.07, 0.05])
    σ_prop  RealDomain(lower = 0,init = 0.03)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @covariates Wt

  @pre begin
    CL = cl * (Wt/70)^0.75 * exp(η[1])
    Vc = tv * (Wt/70) * exp(η[2])
    Ka = ka * exp(η[3])
    Vp = 30.0
    Q  = q
  end

  @dynamics Depots1Central1Periph1

  @derived begin
      cp := @. 1000*(Central / Vc) # We use := because we don't want simobs to store the variable
      dv ~ @. Normal(cp, abs(cp)*σ_prop)
    end
end
PumasModel
  Parameters: cl, tv, ka, q, Ω, σ_prop
  Random effects: η
  Covariates: Wt
  Dynamical variables: Depot, Central, Peripheral
  Derived: dv
  Observed: dv

Some parameters are left free by giving them domains in the @param-block, and one PK parameter (the volume of distribution of the peripheral compartment) is fixed to 20.0.

Simulating the individual observations

The simobs function is used to simulate individual time series. We input the model, the population of Subjects that currently only have dosage regimens and covariates, the parameter vector and the times where we want to simulate. Since we have a proportional error model we avoid observations at time zero to avoid degenerate distributions of the dependent variable. The problem is, that if the concentration is zero the variance in distribution of the explained variable will also be zero. Let's use the default parameters, as set in the @param-block, and simulate the data

param = init_params(mymodel)
Random.seed!(1234)
obs = simobs(mymodel, 
             pop, 
             param, 
             obstimes = 1:1:72)
interactive!(false)
p = sim_plot(mymodel, 
          obs, observations = :dv)
p[1]

Fitting the model

To fit the model, we use the fit function. It requires a model, a population, a named tuple of parameters and a likelihood approximation method.

result = fit(mymodel, 
            Subject.(obs), 
            param, Pumas.FOCEI(), ensemblealg = EnsembleThreads())
Iter     Function value   Gradient norm 
     0     1.002801e+04     1.608331e+02
 * time: 6.198883056640625e-5
     1     1.002794e+04     4.438569e+01
 * time: 0.055644989013671875
     2     1.002778e+04     3.281777e+01
 * time: 0.06789994239807129
     3     1.002723e+04     1.152964e+01
 * time: 0.07953786849975586
     4     1.002646e+04     1.498682e+01
 * time: 0.0889890193939209
     5     1.002549e+04     3.528519e+01
 * time: 0.0987548828125
     6     1.002525e+04     5.396675e+00
 * time: 0.10743594169616699
     7     1.002515e+04     1.900268e+00
 * time: 0.1170189380645752
     8     1.002507e+04     4.798916e+00
 * time: 0.12464499473571777
     9     1.002495e+04     8.321538e+00
 * time: 0.13320684432983398
    10     1.002487e+04     5.994750e+00
 * time: 0.14234185218811035
    11     1.002485e+04     1.719327e+00
 * time: 0.15235590934753418
    12     1.002485e+04     1.420724e-01
 * time: 0.15978598594665527
    13     1.002485e+04     6.651317e-02
 * time: 0.16875195503234863
    14     1.002485e+04     2.056923e-01
 * time: 0.17668485641479492
    15     1.002485e+04     2.617325e-01
 * time: 0.18510890007019043
    16     1.002485e+04     1.458920e-01
 * time: 0.19348788261413574
    17     1.002485e+04     3.011164e-02
 * time: 0.20208191871643066
    18     1.002485e+04     1.971643e-03
 * time: 0.2097618579864502
    19     1.002485e+04     2.117735e-03
 * time: 0.21756601333618164
    20     1.002485e+04     4.095262e-03
 * time: 0.2253720760345459
    21     1.002485e+04     4.289614e-03
 * time: 0.23316693305969238
    22     1.002485e+04     2.136597e-03
 * time: 0.24097490310668945
    23     1.002485e+04     1.319400e-03
 * time: 0.2505300045013428
    24     1.002485e+04     1.183455e-03
 * time: 0.2590649127960205
    25     1.002485e+04     1.183389e-03
 * time: 0.2725410461425781
    26     1.002485e+04     9.994834e-04
 * time: 0.28187084197998047
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                   -10024.853
Number of subjects:                             24
Number of parameters:         Fixed      Optimized
                                  0              8
Observation records:         Active        Missing
    dv:                        1728              0
    Total:                     1728              0

---------------------
           Estimate
---------------------
cl          0.7039
tv          9.774
ka          1.0298
q           0.50035
Ω₁,₁        0.76081
Ω₂,₂        0.048151
Ω₃,₃        0.058265
σ_prop      0.03027
---------------------

Of course, we started the fitting at the true parameters, so let us define our own starting parameters, and fit based on those values

alternative_param = (
    cl = 0.5,
    tv = 9.0,
    ka = 1.3,
    q  = 0.3,
    Ω  = Diagonal([0.18,0.04, 0.03]),
    σ_prop = 0.04)

fit(mymodel, Subject.(obs), alternative_param, Pumas.FOCEI(), ensemblealg = EnsembleThreads())
Iter     Function value   Gradient norm 
     0     2.192254e+04     3.944925e+04
 * time: 4.982948303222656e-5
     1     1.396191e+04     1.736510e+04
 * time: 0.03451800346374512
     2     1.106810e+04     4.365316e+03
 * time: 0.04655599594116211
     3     1.078869e+04     1.310417e+03
 * time: 0.05642580986022949
     4     1.072387e+04     1.278217e+03
 * time: 0.06484198570251465
     5     1.065219e+04     1.840211e+03
 * time: 0.07382488250732422
     6     1.041598e+04     5.547035e+03
 * time: 0.08307480812072754
     7     1.009012e+04     3.261406e+03
 * time: 0.09153485298156738
     8     1.006355e+04     3.884887e+02
 * time: 0.10253190994262695
     9     1.006249e+04     1.069016e+02
 * time: 0.11028385162353516
    10     1.006200e+04     1.048583e+02
 * time: 0.11846089363098145
    11     1.006020e+04     2.220255e+02
 * time: 0.12783384323120117
    12     1.005671e+04     3.359330e+02
 * time: 0.136185884475708
    13     1.005210e+04     3.203693e+02
 * time: 0.1452028751373291
    14     1.004922e+04     1.521796e+02
 * time: 0.1540839672088623
    15     1.004858e+04     3.634242e+01
 * time: 0.16219496726989746
    16     1.004844e+04     3.631365e+01
 * time: 0.1716480255126953
    17     1.004830e+04     4.239529e+01
 * time: 0.17960286140441895
    18     1.004791e+04     9.158692e+01
 * time: 0.187391996383667
    19     1.004701e+04     1.583750e+02
 * time: 0.19507694244384766
    20     1.004501e+04     2.382150e+02
 * time: 0.2052619457244873
    21     1.004146e+04     2.906203e+02
 * time: 0.21512103080749512
    22     1.003672e+04     2.492799e+02
 * time: 0.22508692741394043
    23     1.003337e+04     1.015827e+02
 * time: 0.23441290855407715
    24     1.003270e+04     2.552435e+01
 * time: 0.24299883842468262
    25     1.003265e+04     2.654450e+01
 * time: 0.25171995162963867
    26     1.003259e+04     2.677991e+01
 * time: 0.26029491424560547
    27     1.003242e+04     2.573818e+01
 * time: 0.2692899703979492
    28     1.003209e+04     3.439711e+01
 * time: 0.27742981910705566
    29     1.003152e+04     4.337631e+01
 * time: 0.28592586517333984
    30     1.003090e+04     3.629501e+01
 * time: 0.29424381256103516
    31     1.003058e+04     1.722268e+01
 * time: 0.30318784713745117
    32     1.003049e+04     1.913596e+01
 * time: 0.31111884117126465
    33     1.003046e+04     1.829670e+01
 * time: 0.31977081298828125
    34     1.003041e+04     1.664881e+01
 * time: 0.3277459144592285
    35     1.003029e+04     1.962906e+01
 * time: 0.33602285385131836
    36     1.003001e+04     3.529991e+01
 * time: 0.3442549705505371
    37     1.002941e+04     5.529316e+01
 * time: 0.3540019989013672
    38     1.002836e+04     7.062268e+01
 * time: 0.3624589443206787
    39     1.002697e+04     6.460823e+01
 * time: 0.3717219829559326
    40     1.002578e+04     3.541353e+01
 * time: 0.3831770420074463
    41     1.002521e+04     9.253277e+00
 * time: 0.39461803436279297
    42     1.002508e+04     3.191871e+00
 * time: 0.406674861907959
    43     1.002505e+04     2.254884e+00
 * time: 0.41837596893310547
    44     1.002505e+04     1.944729e+00
 * time: 0.431121826171875
    45     1.002505e+04     1.945441e+00
 * time: 0.44263291358947754
    46     1.002505e+04     1.939199e+00
 * time: 0.45410799980163574
    47     1.002504e+04     1.925990e+00
 * time: 0.46538400650024414
    48     1.002504e+04     1.906824e+00
 * time: 0.4767138957977295
    49     1.002504e+04     1.868390e+00
 * time: 0.4896550178527832
    50     1.002504e+04     1.794271e+00
 * time: 0.501863956451416
    51     1.002503e+04     1.737884e+00
 * time: 0.5245168209075928
    52     1.002501e+04     2.818073e+00
 * time: 0.5375969409942627
    53     1.002496e+04     3.787146e+00
 * time: 0.5500078201293945
    54     1.002490e+04     3.550240e+00
 * time: 0.562675952911377
    55     1.002486e+04     1.667143e+00
 * time: 0.5786800384521484
    56     1.002485e+04     3.962635e-01
 * time: 0.5912308692932129
    57     1.002485e+04     7.931606e-02
 * time: 0.6037249565124512
    58     1.002485e+04     3.453083e-02
 * time: 0.6149590015411377
    59     1.002485e+04     2.767901e-02
 * time: 0.6278948783874512
    60     1.002485e+04     2.767901e-02
 * time: 0.6490998268127441
    61     1.002485e+04     2.767901e-02
 * time: 0.6750028133392334
    62     1.002485e+04     2.769532e-02
 * time: 0.6898720264434814
    63     1.002485e+04     2.769532e-02
 * time: 0.717491865158081
    64     1.002485e+04     1.259742e-01
 * time: 0.7309579849243164
    65     1.002485e+04     6.421964e-02
 * time: 0.7441809177398682
    66     1.002485e+04     3.399018e-02
 * time: 0.7608239650726318
    67     1.002485e+04     1.673434e-01
 * time: 0.7866170406341553
    68     1.002485e+04     3.342486e-01
 * time: 0.8018498420715332
    69     1.002485e+04     4.374753e-01
 * time: 0.8140339851379395
    70     1.002485e+04     3.371089e-01
 * time: 1.5997569561004639
    71     1.002485e+04     1.180698e-01
 * time: 1.6088249683380127
    72     1.002485e+04     1.079942e-02
 * time: 1.6190450191497803
    73     1.002485e+04     1.079942e-02
 * time: 1.6366229057312012
    74     1.002485e+04     1.079942e-02
 * time: 1.656303882598877
    75     1.002485e+04     1.079906e-02
 * time: 1.6758918762207031
    76     1.002485e+04     1.079905e-02
 * time: 1.6911888122558594
    77     1.002485e+04     1.079901e-02
 * time: 1.7079260349273682
    78     1.002485e+04     1.079897e-02
 * time: 1.7238948345184326
    79     1.002485e+04     1.079896e-02
 * time: 1.7403960227966309
    80     1.002485e+04     1.079896e-02
 * time: 1.7574079036712646
    81     1.002485e+04     1.079895e-02
 * time: 1.7761459350585938
    82     1.002485e+04     1.079895e-02
 * time: 1.7963709831237793
    83     1.002485e+04     1.079895e-02
 * time: 1.8134968280792236
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                   -10024.853
Number of subjects:                             24
Number of parameters:         Fixed      Optimized
                                  0              8
Observation records:         Active        Missing
    dv:                        1728              0
    Total:                     1728              0

---------------------
           Estimate
---------------------
cl          0.7039
tv          9.774
ka          1.0298
q           0.50035
Ω₁,₁        0.76079
Ω₂,₂        0.048149
Ω₃,₃        0.058264
σ_prop      0.03027
---------------------

and we see that the estimates are essentially the same up to numerical noise.

To augment the basic information listed when we print the results, we can use infer to provide RSEs and confidence intervals

infer(result)
Asymptotic inference results using sandwich estimator

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                   -10024.853
Number of subjects:                             24
Number of parameters:         Fixed      Optimized
                                  0              8
Observation records:         Active        Missing
    dv:                        1728              0
    Total:                     1728              0

--------------------------------------------------------------------
          Estimate           SE                      95.0% C.I.
--------------------------------------------------------------------
cl         0.7039          0.12541            [0.45809 ;  0.9497  ]
tv         9.774           0.43837            [8.9148  ; 10.633   ]
ka         1.0298          0.052264           [0.92735 ;  1.1322  ]
q          0.50035         0.0011047          [0.49818 ;  0.50251 ]
Ω₁,₁       0.76081         0.2891             [0.19418 ;  1.3274  ]
Ω₂,₂       0.048151        0.013182           [0.022315;  0.073986]
Ω₃,₃       0.058265        0.016079           [0.02675 ;  0.08978 ]
σ_prop     0.03027         0.00051989         [0.029251;  0.031289]
--------------------------------------------------------------------

So as we observed earlier, the parameters look like they have sensible values. The confidence intervals are a bit wide, and especially so for the random effect variability parameters. To see how we can use simulation to better understand the statistical properties of our model, we can simulate a much larger population and try again

pop_big = map(i -> Subject(id = i,
                           events = repeated_dose_regimen,
                           observations = (dv = nothing,),
                           covariates = choose_covariates()),
                           1:100)
obs_big = simobs(mymodel, pop_big, param, obstimes=1:1:72)
result_big = fit(mymodel, Subject.(obs_big), param, Pumas.FOCEI(), ensemblealg = EnsembleThreads())
infer(result_big)
Iter     Function value   Gradient norm 
     0     3.679364e+04     7.630791e+02
 * time: 8.20159912109375e-5
     1     3.679332e+04     1.425848e+02
 * time: 0.16926002502441406
     2     3.679307e+04     1.234952e+02
 * time: 0.21648120880126953
     3     3.679285e+04     5.195340e+01
 * time: 0.7413251399993896
     4     3.679173e+04     7.220420e+01
 * time: 0.7796931266784668
     5     3.679106e+04     4.575628e+01
 * time: 0.8308331966400146
     6     3.679072e+04     1.144429e+01
 * time: 0.869042158126831
     7     3.679068e+04     4.149921e+00
 * time: 0.9115841388702393
     8     3.679067e+04     5.467597e-01
 * time: 0.9500491619110107
     9     3.679067e+04     5.097025e-01
 * time: 0.9821732044219971
    10     3.679066e+04     7.409105e-01
 * time: 1.0126211643218994
    11     3.679066e+04     5.546808e-01
 * time: 1.0436930656433105
    12     3.679066e+04     1.370981e-01
 * time: 1.0739221572875977
    13     3.679066e+04     1.578532e-02
 * time: 1.1063601970672607
    14     3.679066e+04     1.578538e-02
 * time: 1.1588921546936035
    15     3.679066e+04     1.581460e-02
 * time: 1.2266860008239746
    16     3.679066e+04     1.203715e-01
 * time: 1.2714600563049316
    17     3.679066e+04     1.505423e-01
 * time: 1.3377540111541748
Asymptotic inference results using sandwich estimator

Successful minimization:                     false

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                   -36790.662
Number of subjects:                            100
Number of parameters:         Fixed      Optimized
                                  0              8
Observation records:         Active        Missing
    dv:                        7200              0
    Total:                     7200              0

--------------------------------------------------------------------
          Estimate           SE                      95.0% C.I.
--------------------------------------------------------------------
cl         1.215           0.11821            [0.98329 ;  1.4467  ]
tv        10.041           0.27033            [9.5115  ; 10.571   ]
ka         1.0159          0.023039           [0.9707  ;  1.061   ]
q          0.5004          0.00051675         [0.49939 ;  0.50141 ]
Ω₁,₁       0.94724         0.13267            [0.68721 ;  1.2073  ]
Ω₂,₂       0.071929        0.010663           [0.051029;  0.092829]
Ω₃,₃       0.049596        0.006597           [0.036666;  0.062526]
σ_prop     0.029833        0.00024886         [0.029345;  0.030321]
--------------------------------------------------------------------

This time we see similar estimates, but much narrower confidence intervals across the board.

Estimating a misspecified model

To explore some of the diagnostics tools available in Pumas, we can try to set up a model that does not fit out data generating process. This time we propose a one compartent model. The problem with estimating a one compartment model when the data comes from a two compartment model, is that we cannot capture the change in slope on the concentration profile you get with a two compartment model. This means that even if we can capture the model fit someone well on average, we should expect to see systematic trends in the residual diagnostics post estimation.

mymodel_misspec = @model begin
  @param   begin
    cl  RealDomain(lower = 0.0, init = 1.0)
    tv  RealDomain(lower = 0.0, init = 20.0)
    ka  RealDomain(lower = 0.0, init = 1.0)
    Ω   PDiagDomain(init = [0.12, 0.05, 0.08])
    σ_prop  RealDomain(lower = 0, init = 0.03)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    CL = cl * (Wt/70)^0.75 * exp(η[1])
    Vc = tv * (Wt/70) * exp(η[2])
    Ka = ka * exp(η[3])
  end
  @covariates Wt

  @dynamics Depots1Central1

  @derived begin
      cp = @. 1000*(Central / Vc)
      dv ~ @. Normal(cp, abs(cp)*σ_prop)
    end
end
PumasModel
  Parameters: cl, tv, ka, Ω, σ_prop
  Random effects: η
  Covariates: Wt
  Dynamical variables: Depot, Central
  Derived: cp, dv
  Observed: cp, dv
alternative_param_no_q = (
    cl = 0.5,
    tv = 9.0,
    ka = 1.3,
    Ω  = Diagonal([0.18,0.04, 0.03]),
    σ_prop = 0.04)

result_misspec = fit(mymodel_misspec, Subject.(obs), 
                     alternative_param_no_q, 
                     Pumas.FOCEI(), ensemblealg = EnsembleThreads())
Iter     Function value   Gradient norm 
     0     9.617797e+04     1.693060e+05
 * time: 8.0108642578125e-5
     1     2.377470e+04     2.167379e+04
 * time: 0.040180206298828125
     2     2.104349e+04     1.578102e+04
 * time: 0.0628671646118164
     3     1.691859e+04     6.360556e+03
 * time: 0.0836172103881836
     4     1.568972e+04     3.107010e+03
 * time: 0.10488510131835938
     5     1.514405e+04     1.262513e+03
 * time: 0.12670207023620605
     6     1.499042e+04     4.451224e+02
 * time: 0.1480541229248047
     7     1.495861e+04     4.080309e+02
 * time: 0.16955018043518066
     8     1.495423e+04     4.069141e+02
 * time: 0.19259309768676758
     9     1.495225e+04     4.046033e+02
 * time: 0.2149050235748291
    10     1.494680e+04     3.956093e+02
 * time: 0.23790311813354492
    11     1.493489e+04     3.723547e+02
 * time: 0.25963521003723145
    12     1.490880e+04     3.650006e+02
 * time: 0.2830040454864502
    13     1.486507e+04     4.361740e+02
 * time: 0.3068711757659912
    14     1.481223e+04     4.027373e+02
 * time: 0.330794095993042
    15     1.477083e+04     2.426534e+02
 * time: 0.3521871566772461
    16     1.475778e+04     7.890989e+01
 * time: 0.37361812591552734
    17     1.475683e+04     8.094033e+01
 * time: 0.3956911563873291
    18     1.475655e+04     8.196633e+01
 * time: 0.4183080196380615
    19     1.475644e+04     8.206960e+01
 * time: 0.44081997871398926
    20     1.475593e+04     8.192692e+01
 * time: 0.4632279872894287
    21     1.475485e+04     8.077322e+01
 * time: 0.4852120876312256
    22     1.475195e+04     7.638707e+01
 * time: 0.509768009185791
    23     1.474556e+04     6.545525e+01
 * time: 0.5339281558990479
    24     1.473330e+04     7.169812e+01
 * time: 0.5564100742340088
    25     1.471788e+04     6.294130e+01
 * time: 0.5815870761871338
    26     1.470981e+04     2.802554e+01
 * time: 0.6071221828460693
    27     1.470869e+04     1.324312e+01
 * time: 0.6262261867523193
    28     1.470863e+04     1.357644e+01
 * time: 0.6474640369415283
    29     1.470862e+04     1.315681e+01
 * time: 0.6711502075195312
    30     1.470862e+04     1.270641e+01
 * time: 0.6967620849609375
    31     1.470862e+04     1.239475e+01
 * time: 0.7218751907348633
    32     1.470861e+04     1.174411e+01
 * time: 0.7490310668945312
    33     1.470858e+04     1.098416e+01
 * time: 1.4301671981811523
    34     1.470851e+04     1.135417e+01
 * time: 1.4516420364379883
    35     1.470833e+04     1.156782e+01
 * time: 1.4727990627288818
    36     1.470794e+04     1.097780e+01
 * time: 1.4933290481567383
    37     1.470725e+04     8.691388e+00
 * time: 1.5140891075134277
    38     1.470638e+04     6.461460e+00
 * time: 1.5351810455322266
    39     1.470559e+04     4.954152e+00
 * time: 1.556359052658081
    40     1.470520e+04     2.765504e+00
 * time: 1.5781450271606445
    41     1.470508e+04     2.724332e+00
 * time: 1.5994701385498047
    42     1.470504e+04     2.675423e+00
 * time: 1.617218017578125
    43     1.470501e+04     2.616206e+00
 * time: 1.6358060836791992
    44     1.470501e+04     2.563224e+00
 * time: 1.6564819812774658
    45     1.470500e+04     2.535331e+00
 * time: 1.6736981868743896
    46     1.470500e+04     2.527151e+00
 * time: 1.6900320053100586
    47     1.470500e+04     2.521881e+00
 * time: 1.7055530548095703
    48     1.470500e+04     2.512274e+00
 * time: 1.7213730812072754
    49     1.470500e+04     2.498372e+00
 * time: 1.7394030094146729
    50     1.470500e+04     2.473337e+00
 * time: 1.756925106048584
    51     1.470500e+04     2.426461e+00
 * time: 1.7831761837005615
    52     1.470498e+04     2.331011e+00
 * time: 1.8028221130371094
    53     1.470496e+04     2.499398e+00
 * time: 1.8230462074279785
    54     1.470490e+04     2.693471e+00
 * time: 1.8456790447235107
    55     1.470480e+04     2.774137e+00
 * time: 1.8676741123199463
    56     1.470464e+04     2.498396e+00
 * time: 1.8870291709899902
    57     1.470447e+04     1.656331e+00
 * time: 1.908945083618164
    58     1.470439e+04     6.916492e-01
 * time: 1.9297130107879639
    59     1.470438e+04     5.068855e-01
 * time: 1.9506402015686035
    60     1.470438e+04     5.316519e-01
 * time: 1.9687931537628174
    61     1.470438e+04     5.312769e-01
 * time: 1.9852612018585205
    62     1.470438e+04     5.295780e-01
 * time: 2.0021719932556152
    63     1.470438e+04     5.279847e-01
 * time: 2.018195152282715
    64     1.470438e+04     5.243690e-01
 * time: 2.035902976989746
    65     1.470438e+04     5.188250e-01
 * time: 2.0537030696868896
    66     1.470438e+04     5.088026e-01
 * time: 2.0722310543060303
    67     1.470438e+04     4.911665e-01
 * time: 2.0908701419830322
    68     1.470438e+04     7.153115e-01
 * time: 2.110969066619873
    69     1.470437e+04     1.100789e+00
 * time: 2.1326301097869873
    70     1.470436e+04     1.597138e+00
 * time: 2.1561310291290283
    71     1.470434e+04     2.018414e+00
 * time: 2.178997039794922
    72     1.470431e+04     1.914077e+00
 * time: 2.2029471397399902
    73     1.470427e+04     9.915496e-01
 * time: 2.9200031757354736
    74     1.470426e+04     1.919551e-01
 * time: 2.9422202110290527
    75     1.470425e+04     2.509096e-01
 * time: 2.9652931690216064
    76     1.470425e+04     3.312967e-01
 * time: 2.9842851161956787
    77     1.470425e+04     2.571555e-01
 * time: 3.003941059112549
    78     1.470425e+04     8.701710e-02
 * time: 3.025247097015381
    79     1.470425e+04     2.446243e-02
 * time: 3.0460000038146973
    80     1.470425e+04     6.000590e-02
 * time: 3.0664451122283936
    81     1.470425e+04     5.240235e-02
 * time: 3.0837981700897217
    82     1.470425e+04     2.019211e-02
 * time: 3.101785182952881
    83     1.470425e+04     4.286986e-03
 * time: 3.1202361583709717
    84     1.470425e+04     1.288839e-02
 * time: 3.1361382007598877
    85     1.470425e+04     1.114816e-02
 * time: 3.15051007270813
    86     1.470425e+04     4.322042e-03
 * time: 3.1664481163024902
    87     1.470425e+04     9.483345e-04
 * time: 3.1831350326538086
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:               Pumas.FOCE
Log-likelihood value:                   -14704.247
Number of subjects:                             24
Number of parameters:         Fixed      Optimized
                                  0              7
Observation records:         Active        Missing
    dv:                        1728              0
    Total:                     1728              0

---------------------
           Estimate
---------------------
cl          1.0129
tv         21.454
ka     666470.0
Ω₁,₁        0.22529
Ω₂,₂        0.03929
Ω₃,₃        0.053478
σ_prop      0.43937
---------------------

First off, the absorption flow parameters Ka is quite off the charts, so that would be a warning sign off the bat, but let us try to use a tool in the toolbox to asses the goodness of fit. We get these by using the inspect function

res_inspect_2cmp = inspect(result)
res_inspect_1cmp = inspect(result_misspec)
p1 = goodness_of_fit(res_inspect_1cmp)
p2 = goodness_of_fit(res_inspect_2cmp)
1-element Vector{AbstractPlotting.Figure}:
 Figure()
p1[1]
p2[1]

The weighted residuals should be standard normally distributed with throughout the time domain. We see that this is the case for the correctly specified model, but certainly not for the misspecified model. That latter has a very clear pattern in time. This comes from the fact that the one compartment model is not able to capture the change in slope as time progresses, so it can never accurately capture the curves generated by a two compartment model.

Conclusion

This tutorial showed how to use fitting in Pumas.jl based on a simulated data set. There are many more models and simulation experiments to explore. Please try out fit on your own data and model, and reach out if further questions or problems come up.